### This script will create graphs of climate data from multiple sites for multiple months
### Each site month will be produced as an independent graph to aid the visual checking of errors
### Script produced by Collin Storlie on 3rd November 2010

# First define the working directories

in.dir = "/home1/99/jc152199/underlog/"
setwd(in.dir)
out.dir = "/home1/99/jc152199/underlog/maxdiff/"

# Import raw data, check the values, and examine the class of each column of the data frame

air.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv', header=T)
raw.data = read.csv(paste(in.dir,'underlograwdata.csv',sep=''), header=T)
head(raw.data)
str(raw.data)

#Subset data to a few sites only for ease of processing

#raw.sub = subset(raw.data, site=='WU9A2')
#raw.sub2 = subset(raw.data, site=='WU11A2')
#raw.data = rbind(raw.sub, raw.sub2)

#Create columns for year, month, and day, populate them all with NA's

air.data$year = air.data$month = air.data$day = NA

#Populate these new date columns by formatting the already existing complete datetime field first into a class date, then into character (with only the relevant portion of date)
#Then format these characters into numerals for easy sorting during the plotting process

air.data$year = as.numeric(format(as.Date(air.data$date,'%Y-%m-%d'),"%Y"))
air.data$month = as.numeric(format(as.Date(air.data$date,'%Y-%m-%d'),"%m"))
air.data$day = as.numeric(format(as.Date(air.data$date, '%Y-%m-%d'),"%d"))
air.data$airrange = air.data$micro_max - air.data$micro_min

#Aggregate readings based on site, year, month, day and mission using the aggregate function to produce temperature range for under log data


#range.fun = function(x){return(max(x,na.rm=T)-min(x,na.rm=T))}

#tdata.max = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = max)
#tdata.min = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = min)
#tdata.range = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = range.fun)


#Rename columns appropriately and then merge datasets

#names(tdata.max)[6]='underlogmax'
#names(tdata.range)[6]='ulrange'
#names(tdata.min)[6]='underlogmin'

#tdata = merge(tdata.min,tdata.max, by=c('site','year','month','day','mission'))
#tdata = merge(tdata,tdata.range,by=c('site','year','month','day','mission'))

# The following commands will create a smaller dataset aggregating data based on site, day, month, and year BUT NOT MISSION
# Then create plots showing the difference between air max and ul max for each site, month, year, day
daily.max = aggregate(raw.data$logtemp, by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day), FUN=max)
daily.min = aggregate(raw.data$logtemp, by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day), FUN=min)
names(daily.max)[5]='underlogmax'
names(daily.min)[5]='underlogmin'
daily.max$day = as.numeric(daily.max$day)
daily.max$month = as.numeric(daily.max$month)
daily.max$year = as.numeric(daily.max$year)
daily.min$day = as.numeric(daily.min$day)
daily.min$month = as.numeric(daily.min$month)
daily.min$year = as.numeric(daily.min$year) 
tdata = merge(daily.min, daily.max, by=c('site','year','month','day'))
tdata = merge(tdata, air.data, by=c('site','year','month','day'))
tdata$ulrange = tdata$underlogmax - tdata$underlogmin
tdata$maxdiff = tdata$micro_max - tdata$underlogmax
tdata$mindiff = tdata$micro_min - tdata$underlogmin

all.data = tdata

sites = unique(raw.data$site)

#Begin a loop that will plot climate data in chronological order for each site-month, combining all site-months for all years into a single .pdf document

for (xsite in sites) 

{

  sub.tdata = subset(all.data, site==xsite)
  years = unique(sub.tdata$year)
  mission = sort(unique(sub.tdata$mission), decreasing=F) 
  pchs = seq(1:NROW(mission))
  mis.cols = data.frame(mission,pchs,col=c(rainbow(length(pchs))))
  sub.tdata = merge(sub.tdata,mis.cols)
  
  if (nrow(sub.tdata)>0) 
  
  {

  pdf(file=paste(out.dir,xsite,".pdf", sep=""))
  
  cat(xsite,'\n')

    for (xyear in years) 
    
    {
  
      sub2.tdata = subset(sub.tdata, year==xyear)
      #sub2.airdata = subset(sub.airdata, year==xyear)
      
      month.list = sort(unique(sub2.tdata$month), decreasing=F)
      
      if (nrow(sub2.tdata)>0)
      
      cat(xyear,'\n')    
      
      {

        for (xmonth in month.list) 
        
        {
    
          sub3.tdata = subset(sub2.tdata, month==xmonth)

          if  (nrow(sub3.tdata)>1)                                                                                  
          
          {
          
    
          plot(sub3.tdata$day, sub3.tdata$maxdiff, type='p', col=sub3.tdata$col, pch = sub3.tdata$pch, ylim=c(min(sub3.tdata$underlogmax),max(sub3.tdata$underlogmax)+3), xlab = "Day", ylab = "Under Log Max Temp", main=paste(xsite," ",xyear,"-",xmonth,sep=""))
          legend(1,max(sub3.tdata$underlogmax)+2,mis.cols$mission[which(mis.cols$mission %in% sub3.tdata$mission)],pch=mis.cols$pchs[which(mis.cols$mission %in% sub3.tdata$mission)],col=mis.cols$col[which(mis.cols$mission %in% sub3.tdata$mission)])
          
          }
    
        }
    
      } 

    }

  }
                        
  dev.off()

}
